import h5py
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.cluster import DBSCAN
import math
import random
import pylab as P
import matplotlib.patches as mpatches
import scipy 
import numpy.linalg as la
import math as pm
from numpy import genfromtxt
from scipy import stats

PREFIX_U = 'u_r_'
PREFIX_V = 'v_r_'
PREFIX_RHO = 'rho_r_'
PATH = '/home/behollis/DATA/in/ts00050/'
DPATH = '/home/behollis/DATA/out/lockExSt/ts00050/'

DIM = 2
MEM = 950

def readMember(mem):
    ''' Read numpy matrix from disk and return tuple of realization data. '''
    filename_u = PREFIX_U + str(mem) + '.txt'
    filename_v = PREFIX_V + str(mem) + '.txt'
    #filename_rho = PREFIX_RHO + str(mem) + '.txt'
    
    u = genfromtxt(PATH + filename_u)
    u = -1.0*u[0:127,0:127]
    
    v = genfromtxt(PATH + filename_v)
    #rho = genfromtxt(PATH + filename_rho)
    
    return (u, v)# , rho)

def openH5pyFiles():
     #print 'reading hdf5 file...'
    #backward integrated...
    f0 = h5py.File(DPATH+'0.0.hdf5', 'r')
    f1 = h5py.File(DPATH+'0.1.hdf5', 'r')
    f2 = h5py.File(DPATH+'0.2.hdf5', 'r')
    f3 = h5py.File(DPATH+'0.3.hdf5', 'r')
    f4 = h5py.File(DPATH+'0.4.hdf5', 'r')
    f5 = h5py.File(DPATH+'0.5.hdf5', 'r')
    
    #f6 = h5py.File(DPATH+'1.6.hdf5', 'r')
    #f7 = h5py.File(DPATH+'1.7.hdf5', 'r')
    #f8 = h5py.File(DPATH+'1.8.hdf5', 'r')
    #f9 = h5py.File(DPATH+'1.9.hdf5', 'r')
    #f10 = h5py.File(DPATH+'1.10.hdf5', 'r')
    #f11 = h5py.File(DPATH+'1.11.hdf5', 'r')
    
    #f12 = h5py.File(DPATH+'2.12.hdf5', 'r')
    #f13 = h5py.File(DPATH+'2.13.hdf5', 'r')
    #f14 = h5py.File(DPATH+'2.14.hdf5', 'r')
    #f15 = h5py.File(DPATH+'2.15.hdf5', 'r')
    #f16 = h5py.File(DPATH+'2.16.hdf5', 'r')
    #f17 = h5py.File(DPATH+'2.17.hdf5', 'r')
    
    '''
    f0f = h5py.File(DPATH+'0.0.forward.hdf5', 'r')
    #f1f = h5py.File(DPATH+'0.1.forward.hdf5', 'r')
    f2f = h5py.File(DPATH+'0.2.forward.hdf5', 'r')
    #f3f = h5py.File(DPATH+'0.3.forward.hdf5', 'r')
    f4f = h5py.File(DPATH+'0.4.forward.hdf5', 'r')
    #f5f = h5py.File(DPATH+'0.5.forward.hdf5', 'r')
    
    f6f = h5py.File(DPATH+'1.6.forward.hdf5', 'r')
    #f7f = h5py.File(DPATH+'1.7.forward.hdf5', 'r')
    f8f = h5py.File(DPATH+'1.8.forward.hdf5', 'r')
    #f9f = h5py.File(DPATH+'1.9.forward.hdf5', 'r')
    f10f = h5py.File(DPATH+'1.10.forward.hdf5', 'r')
    #f11f = h5py.File(DPATH+'1.11.forward.hdf5', 'r')
    
    f12f = h5py.File(DPATH+'2.12.forward.hdf5', 'r')
    #f13f = h5py.File(DPATH+'2.13.forward.hdf5', 'r')
    f14f = h5py.File(DPATH+'2.14.forward.hdf5', 'r')
    #f15f = h5py.File(DPATH+'2.15.forward.hdf5', 'r')
    f16f = h5py.File(DPATH+'2.16.forward.hdf5', 'r')
    #f17f = h5py.File(DPATH+'2.17.forward.hdf5', 'r')
    '''
    
    print 'finished reading h5py file!'
    
    #iles = ( f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11, f12, f13, f14, f15, f16, f17 )
    #files_f = ( f0f,f1f,f2f,f3f,f4f,f5f,f6f,f7f,f8f,f9f,f10f,f11f, f12f, f13f, f14f, f15f, f16f, f17f )
    
    files = ( f0, f1, f2, f3, f4, f5 )
    #files_f = ( f0f,f1f,f2f,f3f,f4f,f5f,f6f,f7f,f8f,f9f,f10f,f11f, f12f, f13f, f14f, f15f, f16f, f17f )
         
    return files#, files_f

def closeH5pyFiles(files):
    for f in files:
        f.close()
        
        
def getKernelEVFFeatureVector(region_x, region_y, vfield_dict, sx, sy, sz):
    ''' Get region feature vector for velocities in a region. '''
    
    # find plus / minus interval about seed
    xival = int( float(region_x) / 2.0 )
    yival = int( float(region_y) / 2.0 )
    
    DIM = 2
    
    # find (u,v) measurements about seed, that number is the feature length
    range_x = range(sx - xival, sx + xival, 1)
    range_y = range( sy - yival, sy + yival, 1)
    NUM_FEATURES = len(range_x) * len(range_y) * DIM
    
    member_list = list()
    
    all_feat_vecs = list()
    for vf in vfield_dict.keys():
        region_member_vec = list()
        member_list.append(vf)
        for x in range_x:
            for y in range_y:
                u = vfield_dict[vf][0][x, y]
                v = vfield_dict[vf][1][x, y]
                region_member_vec.append( u )
                region_member_vec.append( v )
        all_feat_vecs.append(region_member_vec)
    
    feat_total = np.array( all_feat_vecs )
    MEM = len( vfield_dict.keys() )
    feat_total.reshape( MEM, NUM_FEATURES )
    
   
    
    return feat_total, member_list
 
def readStreamlines(x, y, b_files, f_files, dinteg = 'backward', max_per_file=100):
    ''' loads streamlines thru seed with spatial extents / member extents '''
    
    slines = dict()
    xlst =list()
    ylst = list()
    zlst = list()
    
    cnt = 0
   
    for f in b_files:
        print 'loading file'+str(file)
        
       
            
        groups = f.keys() #use only f for keys, but they should be the same as for ff
        for mem in groups:
            
            if cnt > max_per_file:
                cnt = 0
                break
            
            dir = str(mem).zfill(4) + '/x' + str(x).zfill(3) + '/y' + str(y).zfill(3)
        
            try:
                if dinteg == 'backward':
                    xlst = list(f[dir][0])
                    ylst = list(f[dir][1])
                    zlst = 0.0 * len(list(f[dir][0]))
                '''
                elif dinteg == 'forward':
                    xlst = list(ff[dir][0])
                    ylst = list(ff[dir][1])
                    zlst = 0.0 * len(list(ff[dir][0]))
                else: # 'both'
                    xlst_b = list(f[dir][0][1:])
                    xlst_b.reverse()
                    
                    ylst_b = list(f[dir][1][1:])
                    ylst_b.reverse()
                    
                    zlst_b = [0.0] * len( xlst_b )
                    #zlst_b.reverse()
                    
                    xlst_f = list(ff[dir][0])
                    ylst_f = list(ff[dir][1])
                    zlst_f = [0.0] * len( xlst_f )
                    
                    xlst = xlst_b + xlst_f
                    ylst = ylst_b + ylst_f
                    zlst = zlst_b + zlst_f
                '''
                    
                slines[mem] = (xlst, ylst, zlst)
                
                cnt += 1
                    
            except:
                print 'trying to read MEMBER or SEED that is not in this hdf5 file...'
                
            
   
    return slines

def calcFeatVecCentroid( feats ):
    
    nmems = feats.shape[0]
    nfeats = feats.shape[1]
    sum_feats = np.zeros(shape=(nfeats))
    
    # sum the features
    for idx in range(0, feats.shape[0]):
        sum_feats += feats[idx]
    
    avg_feats = sum_feats / nmems
    return avg_feats

def euclideanDist(pt1, pt2):
    sqdiff = list()
    l1 = len(pt1)
    l2 = len(pt2)
    if l1 is not l2:
        return -1 #not valid operation
    for idx in range(0, l1):
        sqdiff.append( math.pow( (pt1[idx] - pt2[idx]), 2 ) )
        
    total = 0.0
    for val in sqdiff:
        total += val
    
    return math.sqrt(total)
        
if __name__ == '__main__':
    
    files_b = openH5pyFiles()
    MEM = 950
    sx = 29 # seed x coord
    sy = 25 # seed y coord
    
    # gather all of the streamlines
    sl_dict = readStreamlines(sx, sy, files_b, None, dinteg='backward', max_per_file=10)
    
    # gather all of the vector fields
    print len(sl_dict.keys())
    vf_dict = dict()
    for midx in sl_dict.keys():
        print midx
        vfield = readMember( int(midx[3:]) )
        vf_dict[midx] = vfield
        
    rx = 10
    ry = 10
    feat_total, mlist = getKernelEVFFeatureVector(region_x = rx, region_y = ry, 
                                       vfield_dict=vf_dict, sx=sx, sy=sy, sz=0)
    
    '''
    #dbscan on 
    try:
        epsilon = 5.0 #minimum distance between clusters
        db = DBSCAN(eps=epsilon).fit(feat_total)
        
        print 'TOTAL CLUSTERS: ' + str(len(set(db.labels_)))
        
    except:
        print feat_total
        exit()
    
    
    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    core_samples_mask[db.core_sample_indices_] = True
    labels = db.labels_
    
    # Number of clusters in labels, ignoring noise if present.
    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
    
    
    print 'clusters: ' + str( n_clusters_ )
    '''
    
    # find distance of each vfield to each cluster centroid
    mean = calcFeatVecCentroid( feat_total )
    
    sl_dict_dist_to_centroid = dict()
    for idx in range(0, len(mlist)):
        dist = euclideanDist( mean, feat_total[idx] )
        sl_dict_dist_to_centroid[mlist[idx]] = dist
        
    sl_distances = list( sl_dict_dist_to_centroid.values() )
    
    #map sl_distances to 0-1.0 range
    max = max(sl_distances)
    min = min(sl_distances)
    
    for distkey in sl_dict_dist_to_centroid.keys():
        sl_dict_dist_to_centroid[distkey] -= min
        sl_dict_dist_to_centroid[distkey] /= (max - min)
        
    '''
    # find mean distance to each cluster center of cluster members
    
    # find variance (second moment) of streamlines distances to cluster centroid, for each cluster
    
    # render all streamlines in 0-25% range, 26% - 50%, 50% - 75%, 75-100% of gaussian distribution 
    # about each cluster center
    
    colors = list()
    for clus in labels:
        r = random.random()
        g = random.random()
        b = random.random()
        colors.append( (r,g,b) )
    '''
    
    
    fig = plt.figure()
    #ax = fig.gca(projection='3d')
    axes = fig.add_subplot(1, 1, 1, axisbg='blue')
    axes.set_aspect('equal', 'datalim')
    # render streamlines

    
    # find plus / minus interval about seed
    xival = int( float(rx) / 2.0 )
    yival = int( float(ry) / 2.0 )
    
    for slkey in sl_dict.keys():
        sl = sl_dict[slkey]
        midx = mlist.index(slkey)
        #clus = labels[midx]
        #ax.plot(sl[0], sl[1], sl[2], color=colors[int(clus)]) 
        #plt.plot(sl[0], sl[1], color=colors[int(clus)])
        #color_from_map = sl_dict_dist_to_centroid[mlist[idx]]\
        scale = 1.0 - sl_dict_dist_to_centroid[slkey]
        #axes.plot([sx, sy])
        axes.plot([sx-xival, sx+xival ], [sy-yival, sy-yival], color='yellow' )
        axes.plot([sx+xival, sx+xival ], [sy-yival, sy+yival], color='yellow' )
        axes.plot([sx+xival, sx-xival ], [sy+yival, sy+yival], color='yellow' )
        axes.plot([sx-xival, sx-xival ], [sy+yival, sy-yival], color='yellow' )
        axes.plot(sl[0], sl[1], color=(1.0*scale, 1.0*scale, 1.0*scale, 1.0) )
        
         
    #plt.title( 'Label: {0}, Size: {1}'.format(c,size_cluster) )   
    #ax.legend()
    plt.show()
    
    closeH5pyFiles( files_b )
    #closeH5pyFiles( files_f )
    
    print 'finished!'
    
    
    
    
   
    